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The class of materials referred to as ultra high temperature ceramics ( UHTC ), 
including ZrB 2 and Hf B 2 , have high melting points, good chemical stability and 
reasonable mechanical properties. They are candidate materials for applications in- 
volving high temperature environments such as sharp leading edges for hypersonic ve- 
hicles, propulsion systems and refractory crucibles 1-3 . They have also been discussed 
as buffer materials for microelectronics 4 . Unlike most ceramics, these materials are 
distinguished by their high thermal conductivity. High thermal conductivity offers a 
number of advantages for high temperature applications, including improved thermal 
shock resistance and enhanced thermal radiation by efficient distribution of thermal 
energy over available surfaces. 

ZrB -2 and H fB 2 are broadly similar materials, but have nontrivial differences. 
Zr and H f have comparable atomic radii and valence electronic structure, making 
their diboride analogues electronically almost indistinguishable. However, the large 
mass difference between Zr and Hf atoms lead to significant deviations especially 
in the vibrational spectra of the diborides. These differences are expected to affect 
their lattice thermal conductivity as well. The crystal structure of ZrB 2 and HfB 2 
is the AlB 2 -type C 32 with space group P6/mmm. These materials are layered with 
alternating planes of closed packed hexagonal Zr/Hf and open hexagonal, graphitic 
B. Metal atoms are situated directly above/below 6-membered rings in the adjacent 
Boron planes. A detailed ab initio study of the structure and properties of these 
materials was recently presented 6 . 

High thermal conductivity ceramics such as ZrB 2 and H fB 2 have significant con- 
tributions from both electronic and phononic thermal carriers to the total thermal 
conductivity K tof 5 . The electronic component K e can be estimated from the electri- 
cal conductivity a using the Weidemann- Franz (WF) empirical relation, n e = crT, 
where kq is the Lorentz constant (2.45 x 10 -8 1U • H • K -2 ) and T is the temperature. 
The phonon part n p h cannot be measured directly and is usually inferred by subtract- 
ing K e from Ktot- The lattice thermal conductivity can be accessed directly however 
from molecular dynamics (MD) simulations. 

Thermal conductivity measurements for single crystal ZrB 2 at room temperature 
have been reported as 140H 7 /(m ■ K) in the basal direction and 100 W/(m ■ K ) along 
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the c-axis 11 . These single sample results included neither a characterization of defects 
which will reduce K to t nor individual estimates for K e and K p h- Thermal conductivity 
data has not been reported for single crystal HfB 2 although such samples have been 
fabricated 12 . 


Thermal conductivity of poly crystalline ZrB 2 and H fB 2 has been more thor- 
oughly studied 8-10 . Room temperature measurements give 60 W/(m ■ K ) for ZrB 2 
and 104 W/(m ■ K) for HfB 2 . The reduction of K tot relative to single crystals comes 
primarily from the thermal resistance of grain boundaries. Variations of K tot with grain 
size can be substantial 8,10 . For ZrB 2 at 300/1, K e was estimated to be 32 W/(m ■ K) 
using the WF relation, giving 28 W/(m ■ K) for K p h- As a function of temperature, 
K to t f° r polycrystalline ZrB 2 rises gradually up to at least 1800/1 primarily due to a 
monotonically increasing K e . Conversely, K p h decreases as 1/T, saturating to a value 
of 18 W/(m-K) above 1000/1 which is well above the Debye temperature (T D = 750 K) 
for this material. For H fB 2 , Ktot is generally higher than for ZrB 2 although a detailed 
comparison controlling for defect distributions and microstructure has not been done. 
Separate estimates for K p h in HfB 2 have not been made although behavior similar 
to ZrB 2 is expected. Based on these results, we estimate the lattice contribution 
to K tot to be on the order of 40%. This suggests for single crystals that the lattice 
component will be at least 56 W/(m ■ K ) within the planes and 40 W/(m ■ K ) in the 
normal direction. 


Recently, we developed Tersoff style interatomic potentials' for ZrB 2 and HfB 2 . 
In this paper, we report results of atomistic Green-Kubo simulations of the lattice 
thermal conductivity K p h for single crystal ZrB 2 and HfB 2 . Direct measurement of 
K p h is difficult in materials with high K e , therefore, simulations of the lattice compo- 
nent can be useful to isolate that particular thermal conduction mechanism. These 
computations represent the first atomistic simulations for ZrB 2 and HfB 2 and also 
the first nontrivial test of these new potentials. 


3 



I. METHOD 


The Green-Kubo theorem relates the lattice thermal conductivity tensor k %3 to the 
integral of the ensemble average of the time correlation function of the heat current 

Ji(t) 

V r°° 

Kij = Jo W 

where i,j = x,y,z, V is the volume, ks is Boltzmann’s constant and T is the tem- 
perature. The instantaneous heat current can be evaluated from the phase space 
configuration of the system at time t as 

j( (t) = + HCfti • Vj)xij\ (2) 

i i<j 

where e* is the per atom energy of atom i, f l3 is the force on atom i due to neighbor 
j, Vi is the velocity of atom i and x l} = Xi — x 3 where x t is the position vector of atom 

i. 

MD simulations were performed using the LAMMPS package 14 . Systems were 
prepared first with a NPT run to allow the volume to expand, followed by a NVE 
run to equilibrate the system at the desired T . Timsteps ranging from O.lfs to lfs were 
used as required to maintain energy conservation. Larger timesteps were generally 
used at low T and smaller timesteps at higher T. After equilibration, time series 
data for the heat current was collected under NVE. To reduce statistical noise for k, 
long NVE runs were required. For lower temperatures, T = 300/1 for example, runs 
of 10ns with a 1/s timestep were performed. For higher temperatures with smaller 
timesteps, runs with a comparable or larger number of steps were used. The largest 
number of steps considered was 4 x 10' steps. Eight independent simulations were 
performed for each system and at each temperature to control for the randomized 
initial velocities. 

Correlation functions were evaluated from the heat current time series. To increase 
the statistics, correlations were considered in a sliding window through the data set. 
This approach gives errors that increase with the correlation time. The integral of 
the heat current correlation function ( HCCF ) C t3 (t) = (Jj(0) Jj(t)) was evaluated 
using the trapezoid rule. Time correlations up to 200 ps gave well-converged values 
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for the thermal conductivity. 

Recently derived Tersoff style interatomic potentials 7 were used for both ZrB- 2 and 
H fB 2 . For ZrB 2 , two potentials were developed that were denoted “Pot 1” and “Pot 
2” . “Pot 2” however gives lattice thermal conductivity values for k zz that are higher 
than the in-plane component which is not consistent with experiment. “Pot 1” does 
give the correct ordering, and therefore we focus our investigations using it. 

Unit cells contain 12 atoms with 4 Zr/Hf atoms and 8 Boron atoms. Simulation 
cells of size 8 x 8 x 16 were utilized giving 12,288 atoms. Thermal conductivity values 
are well converged with respect to system size at this level. 


II. RESULTS 

Lattice thermal conductivity was obtained for both ZrB 2 and HfB 2 at tempera- 
tures of T = 300/1,600/1, 800/1,1000/1. In Fig 1, the normalized HCCF C xx (t) is 
presented for three different systems: ZrB 2 at T = 300/1 , HfB 2 at T = 300/1 and 
H fB 2 at T — 1000/1 . Analogous results can be obtained for C yy and C zz . A dis- 
tinctive feature of these curves is their rapid and regular oscillation. The oscillation 
frequency and decay for ZrB 2 and H fB 2 at 300/1 is similar, reflecting the broad 
commonalities of these materials. However, for H fB 2 at 1000/1 , there is a signifi- 
cantly stronger damping. These features will be directly reflected in the values for 
k. This oscillatory behavior is different from many HCCF reported in the literature. 
The HCCF for Silicon, for example, decays monotonically without any oscillations 15 . 
Oscillatory HCCF have been observed previously, usually in systems with nontrivial 
mass differences 16 . 

The HCCF oscillations can be analyzed in more detail by examining the Fourier 
transform. In Fig 2, the power spectra for six different HCCF are displayed: the 
C xx and C zz components for ZrB 2 at T = 300/1 and the same components for 
H fB 2 at T = 300/1 and at T = 1000/0 Interestingly, the spectra are peaked 
around single frequencies which can be identified with specific optical phonon modes. 
Phonon spectra for these interatomic potentials have been calculated previously 7 . For 
example, for ZrB 2 at T = 300 K, the power spectra for C xx spikes at to = 364 cm -1 . 
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This frequency corresponds precisely to the in-plane Zr — B optical phonon mode. 
Similarly, the C zz spectra spikes at u = 510cm _1 which matches the out-of-plane 
Zr — B optical phonon mode. Analogous results are shown for HfB 2 at T = 300/1 . 
For HfB 2 at 1000/1 , we see a broadening of the peaks due to anharmonic effects. 

Numerical integration of the HCCF is shown in Fig 3 for ZrB 2 and H fB 2 at 
T = 300/1 as a function of the upper integration limit r. For each material, two 
independent components of the thermal conductivity tensor are shown: K\\ is the 
average of the symmetric in-plane directions k xx and n yy , and k zz is normal to the 
planes. As can be seen, the curves plateau above lOOps indicating largely converged 
integrals. Values for the ZrB 2 lattice thermal conductivity when integrated out to 
200ps are 60 W/(m ■ K) for the in-plane component and 41 W/(m ■ K) in the normal 
direction. For HfB- 2 , the values are 74 W/(m ■ K ) and 63 W/(m ■ K ), respectively. 
For both of these materials, k zz is less than k\\ which agrees with the experimental 
ordering. Values for HfB 2 are also higher than for ZrB 2 . 

Results for the HCCF are relatively free of statistical noise, which is a consequence 
of the very long runs we performed. Other methods have been reported in the liter- 
ature to obtain reliable numbers from shorter runs. For example, fitting the HCCF 
to sums of exponentials has been proposed. However, the correct analytic form is 
not known at this time. Therefore, this approach could introduce errors. Modern 
computer power appears to be adequate to obtain sufficient statistical convergence 
without fitting. 

Next we consider the temperature dependence of the lattice thermal conductiv- 
ity. In Fig 4, we plot both the in-plane acm and the normal k, zz conductivities at 
four different temperatures T = 300 K, 600/1, 800/1, 1000/1 for ZrB 2 and HfB 2 . As 
expected, there is a rapid decrease as a function of T. In this temperature regime, 
anharmonic phonon scattering dominates with the phonon mean free path decreasing 
as 1/T. Fits of the different data sets to 1/T are shown in the graph. All four curves 
converge to similar values for k at high T, which at T — 1000/1 is about 5 W/(m ■ K). 
Compared to experiment however the high T values obtained from our simulations 
for k are probably too low. Lattice thermal conductivity for polycrystalline ZrB 2 for 
example saturates around 18W/(m • K) which can be considered a lower bound for 
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our single crystal results. The low k values at high T most likely reflect a limitation 
in the description of the anharmonic part of the interatomic potentials. 


III. DISCUSSION 

Further insight can be obtained by considering the approximate relation R p h = 
\pCv s \ p h where p is the density, C is the specific heat, v s is the velocity of sound, X p h 
is the phonon mean free path and R p h is the isotropic lattice thermal conductivity. 
Both experimental and computed values for these quantities are given in Table I. The 
specific heat of these materials is comparable with a value of ~ 48 J/ (mol ■ K). The 
velocity of sound can be evaluated using v s = 0.87 \J E / p where E is the bulk modulus, 
giving 5.16/cm/s for ZrB 2 and 3.78 km/ s for HfB 2 based on the experimental values 
while giving 5A5km/s for ZrB 2 and 5.29 km/s for HfB 2 based on values from the 
interatomic potentials. Since v s ~ ZrB 2 is expected to have a higher v s 

than HfB 2 which is reflected in the v s derived from experimental numbers. However, 
the HfB 2 interatomic potential gives a modulus too high by a factor of 2X , thus 
the theoretical values of v s for the two materials are comparable. Estimates for the 
phonon mean free path can be obtained using experimental and simulation results 
for R p h . For ZrB 2 at 300/1, an experimental estimate of 51 Wj(m ■ K) for R p h gives 
\ p h = 11.3 nm while the simulation value of 54 W/ (m ■ K ) for R ph gives X p h = 11.4nm. 
These results are in very good agreement, although an accurate experimental estimate 
for R p h requires further study. For HfB 2 at 300/1 , experimental data is not available. 
The simulation value of 70 W/(m ■ K ) for R p h gives X p h = 14.3nm for that material. 

In this paper, atomistic Green-Knbo simulations for the lattice thermal conduc- 
tivity of ZrB 2 and HfB 2 were performed for a range of temperatures. Interatomic 
potentials for these materials were recently developed and this paper represents the 
first non-trivial test of those potentials. Oscillations in the heat current correlation 
functions can be identified with mixed Zr / H f — B optical phonon modes arising from 
the large metal/Boron mass difference. Agreement with available experimental results 
is very good at room temperature, but is probably too low at higher temperatures 
e.g. T = 1000/1 . 
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p(g/cm?) 

T me lt (K) 

Td(K) 

C(J/ (mol ■ K)) 

E(GPa) 

v s (km/s) 

A P h(nm) 

R ph (W/m ■ K ) 

Zr B 2 Exp. 

6.119 

3518 

750 

48.2 

215 

5.16 

11.3 

51 

MD 





240 

5.45 

11.4 

54 

HfB 2 Exp. 

11.212 

3653 

565 

49.5 

212 

3.78 

N/A 

N/A 

MD 





415 

5.29 

14.3 

70 


TABLE I. Room temperature properties for single crystal ZrB -2 and H fB -2 where p is the 
density, T me i t is the melting temperature, To is the Debye temperature, v s is the velocity 


of sound calculated from v s = 0.87 \jEjp where E is the bulk modulus, X p h is the phonon 
mean free path calculated from R p h = \pCv s \ p h and R p h is the isotropic lattice thermal 
conductivity. 
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FIG. 1 . Normalized heat current correlation function ( HCCF ) for C xx (t) as a function 
of time. Unlike monoatomic systems, the HCCF for ZrB2 and H fB 2 have rapid, regular 
oscillations. 
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FIG. 2. Power spectra of six heat current correlation functions obtained from Fourier 
transformation. Low temperature spectra are peaked around a single frequency identified 
with mixed metal-Boron optical phonon modes. 
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FIG. 3. Integrated HCCF as a function of the upper integration limit r giving both in- 
plane (x-y plane) K|| and out-of-plane (z-direction) k zz conductivities at 300it for ZrF >2 
and H 
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FIG. 4. In-plane km and out-of-plane k zz lattice thermal conductivity as a function of 
temperature for ZrB 2 and H fB- 2 - Curves are fits of the data to 1/T. 
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